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REAL GAS COMPUTATION USING AN ENERGY RELAXATION 
METHOD AND HIGH-ORDER WENO SCHEMES 


PHILIPPE MONTARNAL * AND CHI-WANG SHU t 

Abstract. In this paper, we use a recently developed energy relaxation theory by Coquel and Perthame 
and high order weighted essentially non-oscillatory (WENO) schemes to simulate the Euler equations of 
real gas. The main idea is an energy decomposition into two parts: one part is associated with a simpler 
pressure law and the other part (the nonlinear deviation) is convected with the flow. A relaxation process is 
performed for each time step to ensure that the original pressure law is satisfied. The necessary characteristic 
decomposition for the high order WENO schemes is performed on the characteristic fields based on the first 
part. The algorithm only calls for the original pressure law once per grid point per time step, without the 
need to compute its derivatives or any Riemann solvers. Both one and two dimensional numerical examples 
are shown to illustrate the effectiveness of this approach. 

Key words. Euler equations, real gas, weighted essentially non-oscillatory schemes. 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. In this paper we consider the Euler equations for a real compressible inviscid fluid, 

d t p + div (pu) — 0, t > 0,x e JR d , 
d t pu + div (pu <g> u -f p) = 0, 
d t E + div (( E + p)u) = 0, 

E = \pW \ 2 + pe, 

where the quantities p, u, p, £ and e represent the density, velocity, pressure, total energy and specific internal 
energy, respectively. In addition, there is an equation of state (EOS) of the form p — p(p,e) associated with 
a strictly convex entropy ps(p, e) which satisfies the following entropy inequalities 

(1.2) dtps + div (psu) < 0. 

The pressure law is furthermore assumed to satisfy 

(1.3) p, e (p, € ) > 0* p(p, 0) = 0 and p(p, oo) = oo. 

In the literature research has been done in order to extend classical schemes designed for perfect gas to 
real gases. Collela and Glaz [1] extended the numerical procedure for obtaining the exact Riemann solution 
to a real-gas case, Grossman and Walters [7], Liou, van Leer and Shuen [13] extended the method of flux- 
vector splitting and flux- difference splitting, Montagnc, Yee and Vinokur [16] developed second-order explicit 
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shock-capturing schemes for real gas, Glaistcr [5] presented an extension of approximate linearized Riemann 
solver with different averaged matrices, while Loh and Liou [15] used the generalization of their Lagrangian 
approach (originally proposed for perfect gas) to obtain the real gas Riemann solution. 

Most of the previous proposed methods would require a computation of the pressure law and its deriva- 
tives, or a Riemann solver. This is not only costly but also problematic when there is no analytical expressions 
of the pressure law (for example if we have only table values). 

Recently Coquel and Perthame [2] have introduced an energy relaxation theory for Euler equations of 
real gas. The main idea is to introduce a relaxation of the nonlinear pressure law by considering an energy 
decomposition under the form e = £\ +£2- The internal energy £\ is associated with a simpler pressure law 
Pi (which is taken as the 7-law in this paper), while 62 stands for the nonlinear perturbation and is simply 
convected by the flow. These two energies are also subject to a relaxation process and in the limit of an 
infinite relaxation rate, one recovers the initial pressure law p. 

From this general framework, Coquel and Perthame have also deduced the extension to general pressure 
laws of classical schemes for polytropic gases, which only uses a single call to the pressure law per grid point 
and time step. No derivatives of the pressure law or any Riemann solvers need to be computed. Another 
advantage of their approach is that its implementation does not depend on the particular expression of 
the equation of states. For the first order Godunov scheme, they have shown that this extension satisfies 
stability, entropy and accuracy conditions. Numerical examples have been provided using first order schemes 
by A. In [9], 

For Euler equations of polytropic gas, high order ENO (essentially non-oscillatory) and WENO (weighted 
ENO) schemes [8, 18, 19, 20, 14, 10, 21] have been quite successful in providing high resolution results for 
complicated flow structures and shocks. The idea of ENO schemes is to adaptively choose local stencils so that 
interpolation across a discontinuity is avoided as much as possible. WENO uses a nonlinear combination of all 
stencils to improve upon accuracy and smoothness of numerical fluxes while maintaining the non-oscillatory 
behavior of ENO near a discontinuity. 

The aim of this paper is to study the implementation of this relaxation method with high order WENO 
schemes [10] for real gases. One and two dimensional numerica examples will be given. 

In Section 2 we provide the general framework of the energy relaxation theory of [2] , followed by a short 
description of high order WENO schemes [10]. We then give the details of the construction of the relaxed 
WENO schemes for general gases. In Section 3 numerical examples are given. We start with a description 
of the different equations of states used in this paper, followed by one dimensional shock tube test problems. 
Two dimensional test cases of a smooth vortex, to test the accuracy of the schemes, and of the double Mach 
reflection problem, are then presented. Concluding remarks aie given in Section 4. In the appendices, we 
give the expressions of the Roe matrices for the relaxation system (A) and the two molecular vibrating gas 
(B). 


2. Implementation of the energy relaxation method with WENO. 

2.1. Energy relaxation theory. The principle of the erergy relaxation theory developed by Coquel 
and Perthame [2] is to find a pressure law pi(p A ,e A ) (simpler than p, typically a polytropic law) and an 
internal energy 0(p A ,e A ) so that the system (1.1) and the entrjpy inequality (1.2) can be recovered, in the 
limit of an infinite relaxation rate A (called the equilibrium kmit), from the following system (called the 
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relaxation system): 


( 2 . 1 ) 


d t p x + div ( p x u x ) =0, i>0,x€ H d , 

d t p x u x + div ( p x u x ® u x + p x ) = 0, 

d t E x + div ((E x + p x )u x ) = A p x (e x - <fi(p x ,e x )) , 

dtP X £ 2 + div(p A u A e^) = -A p x ( e £ - <P{p X ,£ 1 )) , 

E \ = ^P X h X \ 2 + P X £ X , 


where pi (p A , e A ) = (71 — 1 )p X £ X with 71 a given constant greater than 1 . One can prove [ 2 ] that the relaxation 
system ( 2 . 1 ) can be supplemented by entropy inequalities under the form 


ftp A E + div(p A Eu A ) < RED a := -Ap A (E, 5 l s lffi - E, £ 2 )(e 2 - <Kp\^)) 

where s\(p,£i) = p ll ~ l /e\ and the specific entropy E denotes an arbitrary function in C^IR^) such that 
pE is convex in (p, p£\, pe 2 ) and that can be written under the form E = E($i(p, £\), £ 2 )- RED A represents 
the Rate of Entropy Dissipation. 

Formally, the original Euler system ( 1 . 1 ) will be recovered at A — ► +00 with 


( 2 . 2 ) 


£ — £l+£2 — £ 1 + 0(p5 £l)i 


provided that we have the following condition (called the consistency condition) 

( 2 - 3 ) P(P,£ 1 + <t>{p,e 1)) = Pi(p,£i) = (71 “ l)pei. 

This last condition can be fulfilled for any given choice of 71 > 1 . 

But in addition to the conservative system ( 1 . 1 ), one also wants to recover at the limit the entropy 
inequality (1.2). The following result, due to Coquel and Perthame [2], gives this last condition under a 
characterization of the admissible 71 . 

Theorem 2.1. Assuming that 71 satisfies 

(2 4 ) 7 i > sup,* r (p, e) where r(p,e) = 1+ 

7 i>sup pe 7 (p,e) where 7 (P> £ ) = pP.c + 7’ 

provided that 71 is finite , we then have 

(i) there exists a (unique) specific entropy E^,^) such that at equilibrium (e = £\ + (j>{p,e\)) 


s(p,e) = S(si(p,£i),^(p,ei)), 

(ii) this entropy is uniformly compatible with the relaxation procedure, i.e .; 

RED X < 0 , for all A > 0. 


2.2. WENO schemes. We use the fifth order WENO scheme in [ 10 ]. For a scalar conservation law 
(2.5) u t +f(u) x = 0 

the derivative f(u) x at the grid point x = Xj is approximated by a conservative flux difference 

(2- 6 ) /(w)i|x=x i ~ ^ (fj+i ~ fj-%) ■ 
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The WENO numerical flux fj+k is computed as follows. For a positive wind direction /'(tc) > 0, we first 
define 3 third order numerical fluxes: 

f - + 1 - 2 ) - 1) f 

(2-7) f] + 1 = -g/(tij_i) + + ^/(Wj+i). 

/ 3 + i = ^/( u j) + ^/( u j+i) - g/(«>+ 2)- 

A third order ENO scheme will result if we choose one of the three third order fluxes in (2.7) adequately, 
according to the size of divided differences [19]. On the other hand, a fifth order linear scheme will result if 
we choose the flux as 


( 2 . 8 ) 




with 

( 2 - 9 ) d! = d 2 = d 3 = — . 

The fifth order WENO scheme results if we choose the numerical flux as 


( 2 . 10 ) 


with Ui defined by 
( 2 . 11 ) 
and 


e; =1 < 


di 


(e + ft) 2 ’ 


01 = Y^(/(w.-2) - 2/(Ui-l) + /(m«)) 2 + 2 ) - 4 /(«i- 1 ) + 3 /(«i)) 2 , 

(2.12) 02 = ^(/(u«-i) - 2 f(ui) + f(u i+ 1)) 2 + i(/(jti_i) - /(w i+ i)) 2 , 

03 = y^(/(tt») - 2/(ui+i) + f{u i+2 )f + 1(3 f(ui) - 4/(u i+ i) + f(u i+2 )) 2 . 

e in (2.11) is taken as 10~ 6 in all our numerical examples, as was done in [10], The weights in (2.11) are 
chosen so that in smooth regions (including at smooth extreira), the WENO flux (2.10) behaves similarly 
as the linear flux (2.8) and is uniformly fifth order accurate. Near shocks, however, the WENO flux (2.10) 
behaves similarly as an ENO flux, in the sense that any stencil crossing a discontinuity has a near zero 
weight. For details of the derivation, see [10, 21]. 

If the wind direction is negative, f'{u) < 0, the procedure is symmetric to the case with f'(u) > 0, with 
respect to the location 2 ^+ 1 - In general, a flux splitting is used 


(2.13) 


f{u) = f+(u) + f (i) 


such that / + (u) has a positive wind direction and / (u) has a negative wind direction: 


(2.14) i-J* W >0, £/-(.)< 0. 

The procedure described above can then be applied to / + (t) and f~{u) separately. The simplest flux 
splitting is the Lax- Friedrichs splitting: 


(2.15) 


/*(«) = \ (/(«)±a«) 
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where 


a = max |/ / (u)|. 

u 

Notice that while first and second order schemes with a Lax-Friedrichs splitting is quite dissipative, higher 
order schemes based on the Lax-Friedrichs fluxes usually give very good results. We use Lax-Friedrichs fluxes 
in this paper. 

For systems of conservation laws in this paper, we use both a component- wise version, where the pro- 
cedure described above is applied to each equation in the system separately, and a characteristic version, 
where locally we apply the procedure described above on characteristic projections. For details about how 
to perform a local characteristic procedure, see for example [18, 19, 10]. 

2.3. Construction of the relaxed WENO scheme. The procedure to solve the Euler system (1.1) 
within the framework of the energy relaxation theory is the following. Given the numerical equilibrium 
solution at the time level t n 

(2.16) p(x,t n ),u(x,t n ),e{x,t n ), 

this approximation is advanced to the next time level t n+1 = t n + At in two steps. 

• First step: relaxation. The two internal energies £i(x,t n ) and £2 (x, £ n ) are obtained by (2.2) and 
the consistency condition (2.3): 

- , ,ns _ pjpjx I * n ),e(s,* n )) 

(2.17) U ’ ’ ( 71 -I )p(x,t") ’ 

e 2 {x,t n ) = £(x,t n ) - Ei(x,t n ). 

Notice that this step involves just one call to the pressure law per grid point and does not involve 
any derivatives of the pressure law or any iterations. 

• Second step: evolution in time. For t n < t < £ n+1 , we solve the Cauchy problem for the relaxation 
system (2.1), with zero on the right side: 

d t p x + div (p A u A ) =0, t > 0, x e IR d , 
d t p x u x + div ( p x u x (8) u x + p x ) = 0, 

(2.18) d t Ei + div ((£^ + Px)u x ) = 0, 

d t p x e 2 + div(p A u A £2) = 0> 

^ A | _ , A |2 i 

El - ~p \u | +/)£„ 

and the initial data 

(2.19) p(x,t n ),u(x,t n ),e 1 (x,t n ),e 2 (x,t n ), 
and we obtain at time t n+1 ~ 

(2.20) p(x, t n+1 ~), u(x, t n+1 “), £i{x, t n+1_ ), £ 2 {x, t n+1 ~). 

At last, we compute the equilibrium solution at time t n+1 by 

p{x,t n+1 ) = p{x,t n+1 ~), 

(2.21) u(x, t n+1 ) = u(x, t n+1 ~), 

e(x,t n+1 ) = £ 1 (x,t n+1 ~) + £ 2 (x,t n+1 ~). 
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Remark 1 . The first step is clearly a relaxation phase, as it is equivalent to the solution of the following 
ODE problem for t >t n 

d t p x = 0, 

. . dtp x u x — 0 , 

(2.22) . ’ . . . . 

d t E x = Xp x {e x -4>{p\e x )), 

dtP X £ 2 = -V A (^2 - <t>(P'\^i)) > 


with initial data at time level t n 


(2.23) p{x,t n ),u{x,t n ),ei(x,t n ).e 2 {x,t n ). 

and to let A — ► +oo. 


We now describe the numerical method we will use for the step of evolution in time. Although our 
numerical results concern both one and two dimensional problems, for simplicity of presentations we shall 
restrict our description to one space dimension. As we are using the finite difference version of WENO 
schemes in [10], extensions to two and more spatial dimensions are simply done dimension by dimension. 
Essentially, the two dimensional code is the one dimensional code with an outside “do loop 1 ’ . 

We have to solve for t n < t < t n+ 1 the following system of four equations 


(2.24) 


d t U + d x F(U) = 0, 

+ initial conditions given by (2.19), 


where 

(2.25) 


U = (p, pu, Ei, pe 2 ) T , 

F(U) = (pu,pu 2 +pi,(Ei + p : )u,pue 2 ) . 


In order to solve the ordinary differential equation 
(2.26) ~U = L(U), 

where L(U) is a discretization of the spatial operator, we use a third-order TVD Runge-Kutta scheme [18]. 

Remark 2. We have two possibilities for the placement of the relaxation step: each Runge-Kutta 
inner stage or each time step . With the first example of section, 3.3, we show that the two approaches give 
nearly identical results in accuracy. Of course the second appr lach is less costly. We thus perform all our 
calculations using the second approach. 

We now discretize the space into uniform intervals of size Ax and denote Xj = j Ax. Various quantities 
at Xj will be identified by the subscript j. 

We use the WENO procedure described in the previous subsection to obtain the spatial operator Lj{U) 
which approximates — d x F{U) at Xj . We have tested several possibilities for the definition of L(U) based 
on WENO schemes. The first one is to use a WENO Lax- Friedrichs scheme with a full characteristic 
decomposition. For this purpose we need to compute a Roe matrix for the system (2.24) and its eigenvalues 
and eigenvectors. The details of this derivation are included in Appendix A. 

The other possibility is to compute the first three components of the numerical flux F 1 i,F 2 i,E? i 

J~r 2 2 J ' 2 

by using a WENO Lax- Friedrichs scheme with a decomposition on the Euler system characteristics and to 
obtain the last numerical flux F 4 k with a scalar WENO Lax Friedrichs scheme. This is possible because 
the first three equations of system (2.24) are independent from the last one. 
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Remark 3. We have also tried to compute the last numerical flux by using a first order scheme specially 
designed in order to preserve the maximum principle for 62 [Hj- But with this approach, we lose the accuracy 
of the high- order WENO scheme also for the other variables . 

Remark 4. In order to make comparisons in the numerical results we have also implemented a WENO 
Lax- Friedrichs scheme with a full characteristic decomposition for a two molecular vibrating gas (see next 
section for a description of the related EOS). For this purpose we need a definition of the corresponding Roe 
average matrix. We give it in Appendix B. For the numerical comparisons for the other real gases we use a 
component- wise WENO Lax- Friedrichs scheme which requires only the computation of the sound velocity 

(2-27) C= V /P ’ P + P ?' 


3. Numerical results. 

3.1. Description of the different equations of states. Wc present here several equations of states 
which we will use in the computation. We find the second one in the paper of In [9], while the third one 
comes from Glaister [4]). 

• Polytropic ideal gas. The equation of states for a polytropic ideal gas (also called perfect gas) is the 
following 

(3.1) p(p,e) = - l)pe. 

Then we have 


( 3 -2) P, p = (7 - 1)£, P,e = (7 - 1 )P- 

Air under normal conditions (p and T moderate enough) can be considered as a perfect gas with 7 = 7/5 — 1.4 
(approximately a mixture of two diatomic molecular species: 20 % of O 2 , 80% of iV 2 ). 

• Two molecular vibrating gas. When the temperature increases the vibrational motion of oxygen and 
nitrogen molecules in air becomes important, and specific heats vary with temperatures. So that one must 
consider the following thermally perfect, calorically imperfect model for two molecular vibrating gas 

(3.3) p(p,e) = rpT(e) 


where the temperature T is given by the implicit expression 

&®vib 


(3.4) 


pe — cJT -f p- 


exp (^ 7 ^) — 1 ’ 

with r = 287.086 J ■ kg C* r = r/( 7 * r — 1), 7 t r = 1.4, G v ib = 10 3 K, a — r. Then we have 

rp 


(3.5) 


P,p = rT(e ), p,, = — 


e'(T(e)) 

• Osborne model. R. K. Osborne from the Los Alamos Scientific Laboratory has developed a quite 
general equation of states in the following form [17] 


(3.6) 


P{p,e) = 


1 

E 4 * <f*o 


(C(oi + a 2C) + E (bo + C(b\ H- & 2 C) + E(cq -f ci£))) 
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where E = p^e and ( = — 1 and the constants po, ai, <* 2 , bo &i, Co, c i> <£o depend on the material in 

question. The typical values for water are p 0 = 10 -2 , a x = 3.84 x 10~ 4 , a 2 = 1.756 x 10 -3 , 6 0 — 1.312 x 1CT 2 , 
bi = 6.265 x 10" 2 , b 2 = 0.2133, c 0 - 0.5132, c x = 0.6761 and = 2 . x 10“ 2 . Then we have 


(3.7) 


“ p 0 (E+<f>o) + 2fl2 0 + E (^1 "" 2 ^C + ^ c l)) j 


E+&; (*>0 + C(^l + ^C) + 2£’(c 0 + CiO) . 


3.2. One dimension cases. 

• ^ramp/e i. ID Riemann problems with perfect gas. 

We consider here two well-known problems which have the following Riemann type initial conditions: 


u(x, 0) 


u l if x < 0 , 

Ur if X > 0. 


The first one is Sod’s problem [ 22 ]. The initial data are 


{pL,u L ,p L ) = (1,0, 1), ( pr, u R , p R ) = (0.125,0,0.1). 


The second one is proposed by Lax [12] with 

(Pl,u l , Pl ) = (0.445,0.698,3.528), {p R ,u Ry p R ) = (0.5,0,0.571). 

Of course for this perfect gas situation there is no need to use the relaxation model in practice. The 
purpose of this test problem is to test the behavior of different relaxation models (different 7 i’s) and differ- 
ent ways of treating the relaxed system (fully characteristic and partially characteristic for the first three 
equations only). 

For this example, a uniform grid of 100 points are used and every 2 points are drawn in the figures. 

We first give, in Table 3 . 1 , a CPU time comparison among the traditional WENO characteristic scheme 
for the perfect gas, and the WENO scheme applied to the relaxation system, both with a fully characteristic 
decomposition and with a partially characteristic decomposition for the first three equations only. The 
calculation is done on a SUN Ultral workstation. We can see that while a fully characteristic decomposition 
is significantly more costly, the partially characteristic decomposition is only slightly more costly than the 
WENO scheme applied to the original perfect gas Euler equations. 


Case 

WENO 

Relaxed WENO with 

Relaxed WENO with 


with characteristic 

full characteristic 

partial characteristic 

Sod Shock 

2.28 

3.49 

2.91 

Lax Shock 

3.32 

4.93 

4.08 


Table 3.1 


CPU time (in seconds) of different schemes for the Sod and Lax shock tube problems for a perfect gas. 


In Figures 3.1 and 3 . 3 , we present the comparison for the Sod’s and Lax’s shock tube problems, of the 
fifth order WENO schemes, applied directly to the perfect gas Euler equations using a characteristic decom- 
position, and applied to the relaxation model with 71 = 3 using only partial characteristic decomposition of 
the first 3 equations. We can see that the results are very close except for the slight over- and under-shoots 
in entropy for the relaxation model calculation. This indicates the feasibility of using the relaxation model. 
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Fig. 3.1. Sod’s shock tube problem with WENO-LF-5 characteristic and relaxed WENO-LF-5 partial characteristic with 
71 = 3.0. 

In Figures 3.2 and 3.4, we present the comparison for the Sod’s and Lax’s shock tube problems, of 
the fifth order WENO schemes. The top left figure compares the full characteristic decomposition for the 
relaxation model, with a partial characteristic decomposition for the first 3 equations only, for 71 = 3. 
We can see that the results are quite close, again indicating the feasibility of using the less costly partial 
characteristic decomposition for the relaxation model. The top right figure compares the effect of different 71 ’s 
in the relaxation model. Apparently bigger 71 corresponds to larger numerical dissipation. This indicates 
that one should always choose the smallest possible 71 subject to stability considerations. The bottom 
figure compares the relaxation WENO results for 71 = 3 and a partial characteristic decomposition, with 
a component-wise WENO scheme applied directly on the original perfect gas Euler equations. Although 
neither uses the correct characteristic information, apparently the relaxation model results are better than 
the component- wise results, especially for the Lax’s problem in Figure 3.4. 

• Example 2. ID Riemann problems with real gases. 

In this example we compute the solutions to the Riemann shock tube problem, for the two molecular 
vibrating gas (3.3)-(3.5) and the Osborne model (3.6)-(3.7), with the following initial conditions in Table 
3.2. 


For this example, a uniform grid of 200 points are used and every 4 points are drawn in the figures. 
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Fig. 3.2. Sod's shock tube problem with WENO-LF- 5 . Comparisons of partial and full characteristic decompositions for 
the relaxation model with 71 = 3 (top left); 71 = 3 and 71 = 30 for the relaxation model with partial characteristic decomposition 
(top right); and the relaxation model with partial characteristic decomposition with 71—3 versus the component- wise WENO 
applied to the original perfect gas Euler equations. 

Also, the “exact solution” in the figures are obtained with the >est scheme using 2000 points. 

We first give a CPU time comparison between the full chara* :teristic decomposition for the original model 
and the partial characteristic decomposition using only the first three equations of the relaxation model, for 
the two molecular vibrating gas model, in Table 3.3. We can see that the partial characteristic decomposition 
for the relaxed model is usually more than twice less costly than the full characteristic version for the original 
system. Although the relaxed model has one more equation, it does not require the computation of the 
complicated derivatives of the EOS. 

In Figure 3.5 we show the comparison of the full characteristic decomposition for the original model and 
the partial characteristic decomposition using only the first three equations of the relaxation model, for the 
two molecular vibrating gas model, with case A initial condition. The results are almost identical, indicating 
that the relaxation model with a partial characteristic decomposition works well with a much reduced cost. 

In Figure 3.6 we show the comparison of the component W!£NO scheme on the original system, and the 
partially characteristic WENO scheme on the relaxed system with 71 = 2.0, for the Osborne gas model with 
case A initial condition. We can sec that the result of the rela ced model is much better, especially for the 
density. This indicates that the relaxation model is a good one for the computation of real gases. 
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Fig. 3.3. Lax’s shock tube problem with WENO-LF-5 characteristic and relaxed WENO-LF-5 partial characteristic with 
71 — 3.0. 

In Figure 3.7 we show the comparison of taking 71 = 10, which satisfies the stability condition (2.4), and 
71 — 2, which satisfies only the second inequality in the stability condition (2.4), for the partial characteristic 
decomposition using only the first three equations of the relaxation model, and the Osborne gas model with 
case A initial condition. We can see that the 71 = 2 results are stable and less dissipative, indicating that in 
practice one does not always have to choose 71 satisfying both inequalities in condition (2.4). 

We have also tested the same problems for the other initial condition cases B, C, D and E. The results 
are mostly similar qualitatively as in case A. To save space we will not present the results here. 

3.3. Two dimensions cases. 

• Example 3 . An isentropic vortex. 

This example is used to verify the accuracy of the relaxation approach, especially the placement of the 
relaxation steps during time stepping. The gas is ideal but we still use the relaxation model. We consider the 
following idealized problem for the Euler equations in 2D: the mean flow is p = 1, p = 1 and (u,u) = (1,1) 
(diagonal flow). We add, to this flow, an isentropic vortex (perturbation in (u, i>) and the temperature 
T = p/p, no perturbation in the entropy S = P/p 1 ): 

(6u, Sv) = L- exp( )(-y,x), ST = ~^~5~ CX P( 1 “ r2 )> ss = 
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Fig. 3.4. Lax’s shock tube problem with WENO-LF-5. Comparisons of partial and full characteristic decompositions for 
the relaxation model with 71 = 3 (top left); 71 = 3 and 71 = 30 for the relaxation model with partial characteristic decomposition 
(top right); and the relaxation model with partial characteristic decomposition with 71 = 3 versus the component- wise WENO 
applied to the original perfect gas Euler equations. 

where (x,y) = (x — 5 , y — 5), r 2 = x 2 + y 2 , and the vortex strength e = 5. See [21]. 

The computational domain is taken as [0, 10, ] x [0, 10], expended periodically in both directions. This 
allows us to perform long time simulation without having to df al with a large domain. 

It is clear that the exact solution of the Euler equation wit h the above initial and boundary conditions 
is just the passive convection of the vortex with the mean velocity. 

In Table 3.4 we show the accuracy result at t = 10 (one time period). We can see that WENO for the 
relaxed model with 71 = 3 gives a somewhat larger error than WENO applied directly to the original system, 
but the order of accuracy is correct. Moreover, to place the relaxation step for each Runge-Kutta inner stage 
or just for each time step seems to give almost identical results. We have thus used the less costly version of 
putting the relaxation step for every time step in all the numerical examples in this paper. 

• Example 4 • Double Mach reflection. 

The computational domain is chosen to be [0, 4] x [0, 1], although only part of it ([0, 3] x [[0, 1]) is shown. 
The reflecting wall lies at the bottom of the computational dor lain starting from x = 1/6. Initially a right- 
moving Mach 10 shock is positioned at (x, y) = (1/6,0) and nakes a 60° angle with the x-axis. For the 
bottom boundary, the exact post-shock condition is imposed tor the part from x = 0 to x = 1/6 and a 
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Case 

State 

Density 

Velocity 

Specific internal 
energy 

A 

Left 

0.066 

0.0 

7.22e6 


Right 

0.030 

0.0 

1.44e6 

B 

Left 

1.40 

0.0 

2 . 22 e 6 


Right 

0.14 

0.0 

2.24e6 

c 

Left 

1.2900 

0.0 

1.95e6 


Right 

0.0129 

0.0 

2.75e6 

D 

Left 

LOO 

0.0 

2 . 00 e 6 


Right 

0.01 

0.0 

2.50e5 

E 

Left 

0.01 

2200.0 

1.44e5 


Right 

0.14 

0.0 

4.00e5 



Table 3.2 



Initial conditions for the test cases 

for real gases 


Case 

WENO 

with characteristic 

Relaxed WENO with 
partial characteristic 

A 

12.68 

5.21 

B 

4.8 

2.63 

c 

12.53 

4.87 

D 

15.0 

5.35 

E 

15.0 

7.84 


Table 3.3 

CPU time (in seconds ) depending on full or partial characteristic decomposition with a two vibrating molecular gas. 


reflective boundary condition is used for the rest. At the top boundary of our computational domain, the 
flow values are set to describe the exact motion of the Mach 10 shock. See [23] for a detailed description of 
this problem. 

First we present the results for a perfect gas. We compare the results using WENO directly on the 
original system [10], and using it on the relaxed model with 71 = 1.5 and 71 — 3.0, in Figure 3.8 for a mesh 
of 480 x 120 points and Figure 3.9 for a mesh of 960 x 240 points. We can see that the relaxed model results 
are quite satisfactory, although a bigger 71 results in some small oscillations. 

Next, we show the results of the same problem with the two vibrating molecular gas. The purpose here 
is to show that the relaxation model based algorithm does work, rather than on the details of the flow with 
more physical models. The results with both a 480 x 120 grid and a 960 x 240 grid are shown in Figure 3.10. 
Comparing with the results in [3], we can see that the main features such as the main shock being closer to 
the bottom boundary, and the shock below the triple point being bent, are also observed here. 

4. Concluding Remarks. We have applied the fifth order WENO schemes to a relaxation model 
to compute the Euler equations of real gases. The algorithm does not depend on the specific form of the 
equation of states and does not need to compute the derivatives of the pressure law. One and two dimensional 
examples are shown to illustrate the accuracy and robustness of the algorithm. 


13 




Fig. 3.5. Case A + two vibrating molecular gas model with WENO-LF-5 characteristic and relaxed WENO-LF-5 partial 
characteristic with 71 = 1.5 


Nb. points 

WENO 

Relaxed WENO 
each time step 

Relaxed WENO 
each R-K step 


LI error 

Accuracy 

LI error 

Accuracy 

LI error 

Accuracy 

20 x 20 

1.07e-2 


1.22e-2 


1.22e-2 


40 x 40 

1.06e-3 

3.3 

2.16e-3 

2.5 

2.17e-3 

2.5 

80 x 80 

6.50e-5 

4.0 

1.77e-4 

3.6 

1.78e-4 

3.6 

160 x 160 

2.09e-6 

4.9 

7.57e-6 

4.6 

7.60e-6 

4.6 


Table 3.4 

LI error and order of accuracy at t - 10 (X period) 


Acknowledgment: We would like to thank Benoit Pertham? for bringing our attention to this problem 
and for helpful discussions. 
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Fig. 3.6. Case A + Osborne gas model with component- wise WENO-LF-5 for the original system and relaxed WENO-LF-5 
partial characteristic with 71 = 2.0. 


Appendix A. Roe matrix for the relaxation system. 

Let us consider two states Ui and U r , then the Roe matrix for the relaxation system ( 2 . 24 ) is the following 


( 


(A.1) 


A(Ui,U . r ) = 


0 1 


0 


(71 ■ 3) T 


u(—Hi + (71 — 1) — ) #i-( 7 i- 1)u 2 


V 


-e 2 u 


7 i“ 

0 


0 \ 


-(71-3 )pu (71-I) 0 

0 

£2 u u ) 


where the averaged state p, u, H\ are defined by 


(A.2) 

P — y/fHy/ ~P t > ^ — otiUi q ' r u rt 

H\ = + a r Hi r , £2 = aj£i ( + o; r £i r 

with 


(A. 3 ) 

= (ei +Pi)/p, 

a . = q = 1 a/ - dr 

1 yfpi+yffr-' Qr 1 Q/ ~ yfpi + y/p^' 


The four eigenvalues of A are 


(A. 4 ) 


a\ — U — C, a 2 = U, 03 = U + C, CI4 — u, 
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Fig. 3.7. Case A + Osborne gas model with the relaxed WENO-LF-5 %:<irtial characteristic with 71 = 10.0 and 71 = 2.0. 


where the averaged sound speed c has the usual expression 


(A.5) 



A set of right eigenvectors can be 


(A.6) 


r\ 


1 1 


1 i 

u — c 


u 

_ 

, r2 = 

u 2 

H — uc 





2 

K ) 


0 / 


1 ^ 

u + c 
H + uc 


r 4 


\ ^ / 


/ 0 \ 
0 
0 


V 1 / 


And the corresponding orthogonal set of left eigenvectors if 



( £(* + §) ^ 


1 — hi \ 


5(61 - f) ) 


^ —£2^1 ^ 


-i(ufe 2 + i) 


b2U 


~\{ub 2- |) 


£ 2 wi>2 

(A.7) /1 = 

62 

> h = 


» h — 

6 2 

» u = 



"2 


— b'2 


2 


-e 2 6 2 


0 / 


V 0 ) 


0 / 


V 1 / 
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Density 


WENO-LF-5 charac. 



Density relaxed WENO-LF-5 partial charac., gammal =1 .5 




Fig. 3.8. Double-Mach reflection, perfect gas, 480 x 120 grid points. 


where 

(A-8) 


(71 - l)u 2 

61 = 2 <? ’ h = 


(71 - 1 ) 
^ 2 


Appendix B. Roe matrix for a two molecular vibrating gas. 

Let us consider two states Ui and t/ r , then the Roe matrix for an Euler system of real gas is the following 
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Density 


WENO-LF-5 charac. 



Density relaxed WENO-LF-5 partial charac., gammal =1 .5 



Density relaxed WENO-LF-5 partial charac., gammal =3.0 



30 contours from 1 .4 to 20.0 Grid:960x240 t=0.I! 

Fig. 3.9. Double-Mach reflection, perfect gas 960 x 240 grid points. 


(see [6] for details) 


(B.l) 


A(U,,U. r) 


0 1 
X + ku 2 /2 (2 - k)u 

\ u(x + «w 2 / 2 — H) H - ku 2 


1 ) 

(1 - I - k)u ) 
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Density Real gas 


relaxed WENO-LF-5 partial charac., gammal =1 .5 



Density Real gas relaxed WENO-LF-5 partial charac., gammal =1 .5 



Fig. 3.10. Double-Mach reflection, two vibrating molecular gas. 


where u, H are the Roe-average values of the velocity and the total specific enthalpy (H = e + 1/2 u 2 + p/p) 


(B.2) 

(B.3) 


- yfPl'U'l T y/Pr^r 

U ~ V^+V^ ’ 

_ y/pAh + s/frHr 
y/Pl+y/P^ 


and x and k are two parameters which must satisfy 


(B-4) 


A p = kA pe + \Ap 


with A p — p r - pi , Ape = p T e r - pie t and A p = p(p r , e r ) - p(pi,ei). 


The definitions for k and \ proposed by In [9] for a two molecular vibrating gas are 


(B-5) 


(B.6) 



r(r(e,)-r(e,)) 

e r-ej 

l/P,«(Pl.e) i P,«(Pr,e) 
2 ^ PI P r 




if e r / ei, 
if e r —e i— e, 


< 


K= < 


Ap—nApe 
A p 


sC pAp< £ i ) - ^pAp^i) + pAp^t) - ^pAp^t)) 
= \r{T{ei) - e,r(e t ) + T(e r ) - e r T'(e r )) 


if p r ± pi, 


if p r = Pi = p. 
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The definitions of the eigenvalues and right and left eigenvectors are easy to obtain and are omitted 
here. 
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